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Abstract The characterization of plasticity, robustness, and evolvability, an 
important issue in biology, is studied in terms of phenotypic fluctuations. By 
numerically evolving gene regulatory networks, the proportionality between 
the phenotypic variances of epigenetic and genetic origins is confirmed. The 
former is given by the variance of the phenotypic fluctuation due to noise 
in the developmental process; and the latter, by the variance of the pheno- 
typic fluctuation due to genetic mutation. The relationship suggests a link 
between robustness to noise and to mutation, since robustness can be defined 
by the sharpness of the distribution of the phenotype. Next, the proportion- 
ality between the variances is demonstrated to also hold over expressions 
of different genes (phenotypic traits) when the system acquires robustness 
through the evolution. Then, evolution under environmental variation is nu- 
merically investigated and it is found that both the adaptability to a novel 
environment and the robustness are made compatible when a certain degree 
of phenotypic fluctuations exists due to noise. The highest adaptability is 
achieved at a certain noise level at which the gene expression dynamics are 
near the critical state to lose the robustness. Based on our results, we revisit 
Waddington's canalization and genetic assimilation with regard to the two 
types of phenotypic fluctuations. 

Keywords Robustness, Fluctuation-Response Relationship, Evolution, 
Genetic Variance 



1 Introduction 

In evolutionary biology, plasticity and robustness are considered the basic 
characteristics of phenotypes; these characteristics have been widely dis- 
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cussed for decades. In general, phenotypes are shaped from genotypes as 
a result of developmental dynamicfO], under an environmental condition. Al- 
though these dynamics are determined by genes, they can be stochastic in 
nature, owing to some noise in the developmental process; further, these dy- 
namics also depend on the environmental condition. 

Plasticity refers to the changeability of the phenotype against the en- 
vironmental change. Through developmental dynamics, the influence of the 
environment is amplified or reduced[l-6]. In other words, plasticity concerns 
how the developmental dynamics are affected by the environmental change. 

Of course, the phenotype depends on the genotype. This changeability 
against genetic change is called evolvability and is related to the sensitivity of 
developmental dynamics against genetic change. In this sense, both plasticity 
and evolvability represent the responsiveness against external perturbations 
and, in simple terms, a sort of "susceptibility" in statistical physics. 

Robustness, on the other hand, is defined as the ability to function against 
possible disturbances in the system [7-13]. Such disturbances have two dis- 
tinct origins: non-genetic and genetic. The former concerns the robustness 
against the stochasticity that can arise during the developmental process, 
while the latter concerns the structural robustness of the phenotype, i.e., its 
rigidity against the genetic changes produced by mutations. If the variance 
of phenotype owing to disturbances such as noise in gene expression dynam- 
ics is smaller, then the robustness is increased. In this sense, the variance 
of phenotypic fiuctuations serves as a measure of robustness, as has already 
been discussed [111. 

Now, there exist certain basic questions associated with robustness and 
evolution[l-8]. Does robustness increase or decrease through evolution? If it 
increases, the rigidity of the phenotype against perturbations also increases. 
Consequently, the plasticity, as well as evolvability, may decrease with evolu- 
tion. If that is the case, then the question that arises is: How can the plasticity 
needed to cope with a novel environment be sustained? 

The decrease in plasticity and increase in robustness with evolution was 
actually observed in laboratory experiments under fixed environmental and 
fitness conditions and was also confirmed through numerical experiments. 
In a simulated evolution of catalytic reaction and gene regulatory networks 
under a given single fitness condition, the fluctuation decreases through the 
course of evolution |14[[T5] . Such networks evolve to reduce the fluctuation by 
noise, for a sufficient noise level. Through this evolution, robustness against 
noise increases, leading to the decrease in phenotypic plasticity. However, as 
a system is more and more adapted to one environment, the phenotype fitted 
for it would lose the potentiality to adapt to a novel environmental condition. 

In the present organisms, however, neither the evolutionary potential nor 
the phenotypic fluctuation vanishes. Even after evolution, the phenotype in 
question is not necessarily fixed at its optimal value, but its variance often 
remains rather large. There can be several sources for the deviation from the 

^ Here, "development" refers to a dynamic process that shapes the phenotype 
and is not restricted to multicelhilar organisms. Developmental dynamics are also 
observed in unicellular organisms, for instance, gene (protein) expression dynamics 
by regulatory networks. 
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optimal state, which is neglected in the idealized numerical and laboratory 
experiments with a single fitness condition. One of the most typical factors 
for this deviation is environmental variation. With environmental change, the 
phenotype for producing the fittest state is not fixed but may vary over gen- 
erations. Then, we are faced with the following questions: First, under envi- 
ronmental variation, can robustness and plasticity coexist through the course 
of evolution? Second, are the phenotypic fluctuations sustained, to maintain 
the adaptability to environmental changes? Third, how are the opposing fea- 
tures of robustness and adaptability to the new environment compromised? 
Finally, is there a noise level optimal for achieving both adaptability and 
robustness? We address these questions in the present paper. 

2 Background: Isogenic phenotypic fluctuation and evolution 

speed 

As for the fluctuation, there is an established relationship between the evolu- 
tionary speed and the variance of the fitness, namely, the so-called fundamen- 
tal theorem of natural selection proposed by Fisher [T^[^[^ : this theorem 
states that the evolution speed is proportional to Vg, the variance of fit- 
ness due to genetic variation. Besides the fitness itself, any phenotypes are 
generally changed by the genetic variation. Now, relationship between the 
evolution speed of each phenotype with its variance due to genetic change is 
established as Price equation, formulated through the covariance between a 
phenotype and fitness i22u23j . Statistical-mechanical interpretation of Fisher's 
theorem in terms of the fluctuation-response relationship |24| or fluctuation 
theorem[3S] was also discussed. 

Besides the genetic change, there are sources of phenotypic fluctuations 
even among isogenic individuals, and under a flxed environment. In fact, 
recent observations in cell biology show that there are relatively large fluctu- 
ations even among isogenic individuals. The protein abundances over isogenic 
individual cells exhibit a rather broad distribution, i.e., the concentration of 
molecules exhibits quite a large variance over isogenic cells 27, 28, 29, . This 
variance is due to stochastic gene expressions and to other external pertur- 
bations. Furthermore some of such phenotypic fluctuations are indeed tightly 
correlated with the fltness, and hence the fltness also shows (relatively large) 
isogenic fluctuations. For example, large fluctuations in the growth speed (or 
division time) are observed among isogenic bacterial cells |30l[5T| . 

Hence, there are two sources of variations for the fitness or phenotype: ge- 
netic and epigenetic. The former concerns structural change in developmental 
dynamics with genetic change, and the latter concerns the noise during the 
gene expression dynamics. As mentioned, the relationship between the for- 
mer variance with the evolution speed was established as Fisher's theorem 
(for fitness) and Price equation (for general phenotype). Then, is there any 
relationship between the latter variance with the evolution speed? 

At a glance, the latter variance might not seem to be relevant to evo- 
lution, since this epigenetic change itself due to noise is not inherited to 
the offspring, in general. This is not the case. Here, it should be noted that 
the degree of variance itself is a nature of developmental dynamics to shape 
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the phenotype, and thus depends on genotype, and can be inherited. Hence, 
there may exist some correlation between evolution speed and this isogenic 
phenotypic variance, which is denoted as Vip hereQ 

Indeed, from evolution experiments for increasing the fluorescence in an 
inserted protein in bacteria|17j and also from numerical experiments for 
evolving reaction networks or gene regulatory networks [1511 14| to increase 
a given fitness, we have observed that 

Vip oc Evolution Speed, 
where the evolution speed is defined as the increase of the fitness. In the 
experiment, for each mutant bacteria, average fluorescence of isogenic cells 
was measured and that with highest average fluorescence was selected. Si- 
multaneously, the variance of isogenic bacterial population was measured 
to obtain Vip. The same procedure was applied to measure the evolution 
speed and the fitness variance in numerical evolution to increase a given 
fitness. Interestingly, these experiments support the above relationship, at 
least approximately. Some 'phenomenological explanation' was proposed by 
assuming a Gaussian-type distribution P{x] a) of the fitness (phenotype) x 
as parameterized by a "genetic" parameter, a. and a linear change in the 
peak position of x against a |I6[[T71[T5] . It should be noted, however, that the 
argument based on the distribution a) is not a "derivation" but rather 
a phenomenological description. Indeed, the description by P{x] a) itself is 
an assumption: for example, whether a genotype is represented by a scalar 
parameter a is an assumption. 

Now, considering the established Fisher's relationship, the above relation- 
ship suggests the proportionality between Vip and Vg through an evolutionary 
course. Indeed, this relationship was confirmed from several simulations of 
models. Again, this relationship is not derived from established relationships 
in population genetics. Indeed, the proportionality between the two is not 
observed in the first few generations, but observed after robust evolution 
preserving a single peak in the fitness is progressed. Considering this obser- 
vation, we previously discussed the relationship, by postulating evolutionary 
stability of the distribution over phenotype and genotype 14,15,35 . 

The relationship between Vip and Vg also suggests a possible link between 
developmental robustness to noise and evolutionary robustness to genetic 
changes (mutation). For this link, we first note that the two types of variances 
Vg and Vip lead to two kinds of robustness: rigidity of the fitness (phenotype) 
against genetic changes produced by mutations and robustness against the 
stochasticity that can arise in an environment or during the developmental 
process. When the fitness (phenotype) is robust to noise in developmental 
process, it is rather insensitive to the noise, and therefore, its distribution 
is sharper. Hence, the (inverse of the) variance of isogenic distribution, Vip, 
gives an index for robustness to noise in developmental dvnamics 1 1 4111 5) . On 
the other hand, if Vg is smaller, the phenotype is rather insensitive to ge- 

^ This is not standard terminology. Phenotypic variance by non-genetic origins 
is often termed as environmental variation 14- Here, however, the source of the 
variance is not necessarily environmental change but, primarily, noise in the devel- 
opmental process; hence, we use a different term. Note th at the total phenotypic 
variance is Vip + Vg under a certain ideal condition [20ll21j . if they are added inde- 
pendently. 
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netic changes, implying higher genetic (or mutational) robustness. Hence the 
correlation between Vip and Vg implies correlation between the two types 
of robustness. Indeed, previous simulations suggest that robustness to noise 
fosters robustness to mutation. Congruence between evolutionary and de- 
velopmental robustness was also discussed by Ancel and Fontana for the 
evolution of RNA^. 

To close this section, a brief remark on the measurement of Vg is given 
here. Because the fitness (phenotype) distribution exists even in isogenic 
individuals, the variance of its distribution over a heterogenic population 
includes both the variance among isogenic individuals and that due to genetic 
variation. To distinguish between the two contributions, we first measure the 
average fitness (phenotype) over isogenic individuals and then compute the 
variance of this average over a heterogenic population. This variance will only 
be attributed to genetic heterogeneity. This variance is denoted by V^. 



3 Our Standpoint and Model 

3.1 Dynamical-systems approach to evolution-development relationship 

To discuss the evolution, the fitness landscape as a function of genotype is 
often adopted following Wright's picturei32j. Energy-like fitness function is 
assigned to a genotype space, i.e.. Fitness = f (Genotype). Here the evolu- 
tion is discussed as a hill-climbing process through random change in geno- 
type space and selection. Although this viewpoint has been important in 
evolutionary studies, another facet in evolution has to be also considered, 
to discuss the phenotypic plasticity and evolution-development relationship, 
that is genotype-phenotype mapping shaped by the developmental dynamics. 
Here we focus on this evolution-development relationship. 

Note that the fitness is a function of (some) phenotypes (say a set of 
protein abundances). The phenotypes are determined by developmental dy- 
namics (say gene expression dynamics), whose rule (e.g., set of equations or 
parameters therein) is governed by genes. Thus this evolution-development 
scheme is represented as follows: 

(i) Fitness= F(phenotype) 

(ii) Phenotype determined as a result of developmental dynamics 

(iii) Rule of developmental dynamics given by genotype 
According to (ii) and (iii), genotype-phenotype mapping is shaped, which 

is not necessarily deterministic. As discussed, the phenotypes from isogenic 
distribution are distributed, since the dynamics (ii) involve stochasticity (in 
gene expression). Indeed, with this stochasticty, reached attractors by the 
dynamics are not necessarily identical, and accordingly the fitness may vary 
even among the isogenic individuals, which lead to the variance Vip. If the 
dynamics to shape the phenotype is very stable, the variance Vip is smaller. 

^ According to the conventional terminology in population genetics, this variance 
is referred to as "additive" genetic variance; The term "additive" is included, to 
remove the variance due to sexual recombination. Here we do not discuss the influ- 
ence of recombination and, thus, do not need to distinguish genetic variance from 
the additive one. 
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Further, the degree of the change in phenotype by the genetic change depends 
both on (ii) and (iii), and Vg depends on the sensitivity of the phenotype to 
the change in the rule of the dynamics. With the evolutionary process the 
genotype, i.e., the rule of the dynamics, changes, so that Vip and Vg change 
with the evolution. 

The developmental dynamics, in general, involve a large number of vari- 
ables (e.g., proteins expressed by genes), and are complex. Accordingly, genotype- 
phenotype mapping is generally complex, and the mapping from genotype 
to the fitness is not simple, even if the function in (i) is simple. In several 
studies with adaptive fitness landscape, complexity in the fitness landscape 
has been taken into account, while the developmental dynamics (ii) are not. 
We take a simple fitness function (say, the number of expressed genes), but 
instead take complex developmental dynamics into account, to discuss the 
phenotypic plasticity in developmental and evolutionary dynamics. 

If we are concerned only with the fitness landscape (i), we can introduce 
a single 'energy'-like function for it and formulate the adaptive evolution in 
terms of standard statistical mechanics. Here, however, we need to consider 
the dynamics to shape the phenotype (ii) [51 [T51[T^IT5] . If we adopt statistical- 
mechanical formulation, we need two 'energy'-like functions, one for fitness 
and the other for Hamiltonian for the development dynamics, as is formulated 
by two-temperature statistical phvsics [331134) . 

So far, such study on evolution-development (so called 'evo-devo') lacks 
mathematically sophisticated formulation as compared with the celebrated 
population genetics, and thus we sometimes have to resort to heuristic studies 
based on numerical evolution experiments, from which we extract 'generic' 
properties. We then provide some plausible arguments (or phenomenologi- 
cal theory), while mathematical or statistical-physical formulation has to be 
pursued in future. 



3.2 Gene expression network model 

Following the argument in the lase section and to discuss the issue of plas- 
ticity and robustness in genotype-phenotype mapping, we adopted a simple 
model[14','39''40' for gene expression dynamics with a sigmoid input-output 
behavior[36,37,38 ; it should be noted, however, that several other simula- 
tions in the form of other biological networks will give essentially the same 
result. In this model, the dynamics of a given gene expression level, Xi, is 
described by the following: 



where Jij = —1,1,0. The noise term rii{t) is due to stochasticity in gene 
expression. For simplicity, it is taken to be Gaussian white noise given by 
< rii{t)rij{t') >= 5ijS{t — t'), while the (qualitative) results to be discussed 



M 
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below is independent of the choice 0- The amphtude of noise strength is given 
by cr, which determines stochasticity in gene expression. M denotes the total 
number of genes; and k, the number of target genes that determine fitness. 
Here, the function f{x) represents a threshold function for gene expression. 
Previously, we adopted the model (1) with 

/(x) -tanh(/3a;), (2) 

where the gene expression is "on" if x > and "off' if a; < 0. 

To discuss the environmental condition, however, it is relevant to intro- 
duce an "input" term for expressing some genes. To make this input effective, 
we modify the function (3) as 

f{x) = l/{exp(~(3(x - e,) + l) + 5 (3) 

Here, x is positive and is scaled so that the maximal level expression is 
~ 1; (5 takes a small positive value corresponding to a spontaneous expres- 
sion level. If the input sum from other genes Yl^/ Jij^j to g^ne i exceeds the 
threshold 9i, the gene (protein) i is expressed so that ~ 1, as / approaches 
a step function for sufficiently large f3{> 1), which corresponds to the Hill co- 
efficient in cell biology. Now, if all initially smaller than the threshold 

they remains so li 5 < 9i. Hence, to generate some expression patterns, 
some input term is needed. Here, we introduce "input" genes j = 1, ...kinp, 
in which Xj is set at 

Xj ^ Ij >0 {j = l,2,...,hnp), (4) 

where Ij is an external input to a set of "input genes." These inputs are 
needed to express genes; otherwise, the expression levels remain sub-threshold. 
In this case with eq. (3), the gene is expressed if Xi > 9i, The initial condition 
is given by a state where none of the genes are expressed, i.e., Xi ~ 0. 

Next, we set the fitness for evolution. The fitness, F, is determined by 
whether the expressions of the given " target" genes are expressed after a suf- 
ficient time. This fitness condition is given such that k target genes are "on" 
(expressed), i.e., Xi > 6i for M — k < i < M. Because the model includes a 
noise component, the fitness can fiuctuate at each run, which leads to a distri- 
bution in F and Xi , even among individuals sharing the same gene regulatory 
network. For each network, we compute the average fitness F over L runs as 
well as the variance. This variance over L identical individuals (having the 
identical networks Jij) due to noise is Vip, as it represents the fiuctuation of 
the fitness over isogenic individuals (isogenic phenotypic variance). Vg, to be 
discussed below, is the variance of F over N heterogenic individuals (having 
different networks Jij). 

* Indeed, multiplicative noise depending on Xi, as well as a stochastic reaction 
model si mu lated with the use of Gillespie algorithm gives a qualitatively same 
behavioral]. Of course the magnitude of the phenotype variance as well as the 
threshold noise level for error catastrophe, to be discussed below, depends on the 
form of noise. However, relationship of such threshold noise level with Vg/Vip, as 
well as the proportionality between Vip and Vg does not depend on the specific 
form of noise. 
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Here Vip is not the variance over time in a single run, but the variance 
over runs for identical individuals. With developmental process by gene ex- 
pression dynamics from a given initial condition, the gene expression pattern 
{x{j);j = 1,..,M) reaches an attractor. Reached attractors can be differ- 
ent by each run, due to the stochasticity in expression dynamics during the 
transient time steps to reach the attractor. The fitness is computed on this 
attractor, which is also distributed as a result of noise. Once the attractor 
is reached, the fluctuation of the fitness around it is negligible, since the fit- 
ness is not given directly by Xi but defined through a threshold function of 
Xi — 9i. The variance is smaller as more runs result in the same attractor (or 
attractors of the same fitness) under noise. If gene expression dynamics with 
such global attraction to the same-fitness attractor are shaped through the 
evolution, the variance Vip is smaller. 

Now, at each generation, there are N individuals with slightly different 
gene regulatory networks {Jij}, and their average fitness F may differ by 
each. Among the networks, we select the ones with higher fitness values. 
From each of the selected networks, Jij is "mutated," i.e., Jij for a certain 
pair i, j selected randomly with a certain fraction switches to one of the values 
among ±1,0. Ns{< N) networks with higher F values are selected, each of 
which produces N/Ng mutants. We repeat this selection-mutation process 
over generations. (For example, we choose N = L = 500 or 700 for most 
simulations, and Ng — N/A; the conclusion, to be shown below, does not 
change as long as these values are sufficiently large. We use /3 = 7, 7 = .1, 
M = 64, and k = 8 and initially choose Jij randomly. The probability is 
taken to be 1/4 for Jy = ±1 and 1/2 for for most simulations below, but 
the results below are independent of this choice [ID]. 

Now, we have two types of variances. Besides Vip, Vg is defined as the 
variance of F over the N individuals having different genes (gene regulatory 
networks). This shows the variance due to genetic change. As Vg is decreased, 
the fitness becomes insensitive to the genetic change, i.e., the robustness to 
mutation is increased. On the other hand, as Vip is decreased, the robustness 
to noise is increased, since Vip measures the variance due to noise in a de- 
velopmental process. (It should be remarked that the absolute value of the 
fitness F is not important, since the model behavior is invariant against the 
transformation F c x F. Accordingly the absolute values of Vip and Vg 
are not important, but the ratio between the two or relative change of the 
variances over generations is important. On the other hand, each expression 
x{i) is already scaled so that the maximal expression is ^ 1). 

4 Confirmation of Relationsiiips in Plienotypic Fluctuations 

4.1 Proportionality between the phenotypic variances by noise and by 
mutation 

From the numerical simulation of the model, we have confirmed the following 
results. 

(1) There is a certain threshold noise level Cc beyond which the evolution 
of robustness progresses, so that both Vg and Vip decrease. Here, most of the 
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Generation 

Fig. 1 Evolutionary time course of the average fitness < F >. The average of 
mean fitness < F > over all A'^ = 500 individuals that have different genotypes 
(i.e., networks Jij) in each generation is plotted. The mean fitness of each genotype 
is computed from L = 500 runs. The plotted points are for different values of 
noise strength, a =0.005, 0.03, 0.06, with different colors. For Figs. 1-3, we choose 
M = 64, fc = 8, and kinp = 8 with Ij = 1, while Oi is distributed in [0.1, .0.3]. 



individuals take the highest fitness value. In contrast, for a lower noise level 
(7 < CTc, mutants that have very low fitness values always remain. Several 
individuals take the highest fitness value, whereas the fraction of individuals 
with much lower fitness values does not decrease; hence, Vg remains large 
(see Figs.l and 2). 

(2) At around the threshold noise level, Vg approaches Vip. For a < ac, 
Vg ^ Vip holds, whereas for a > (Tc, Vip > Vg is satisfied. For robust evolution 
to progress, this inequality is satisfied. 

(3) When the noise is larger than this threshold, the two variances de- 
crease, while Vg oc Vip is maintained through the evolution course. Hence, 
the proportionality between the two variances is confirmed. 

Why does the system not maintain the highest fitness state under a small 
phenotypic noise level with a < del Indeed, the dynamics of the top- fitness 
networks that evolved under such low noise levels have distinguishable fea- 
tures from those that evolved under high noise levels. It was found that 
for networks evolved under cr > CTc, a large portion of the initial conditions 
reached attractors that give the highest fitness values, whereas for networks 
evolved under cr < dc, only a tiny fraction (i.e., in the vicinity of the all-off 
states) reached such attractors. 

In other words, for a > ac, the "developmental" dynamics that give a 
functional phenotype have a global, smooth attraction to the target. In fact, 
such types of developmental dynamics with global attraction are known to be 
ubiquitous in protein folding dynamics [42ti44j . gene expression dynamics [43], 
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Fig. 2 Relationship between Vg and Vip, the variances of the fitness. Vg is com- 
puted from P{F), the distribution of mean fitness F and Vip is computed from the 
isogenic variance of the fitness over L — 500 runs, for each genotype, and then these 
values are averaged over all existing individuals. (We also confirmed the overall re- 
lationship by using the variance for a gene regulatory network that gives the peak 
fitness value in P{F)). The plotted points are over 140 generations, a = 0.005 (□), 
0.03 (*), 0.05 (x) and 0.06 {+). For a > Uc « 0.02, both the variances decrease 
with generatons, so that the right-top is the first generation and left-bottom is the 
140th generation. After initial few generations both decrease roughly in proportion, 
while a.t a — 0.03, the deviation from the proportionality is larger possibly because 
the noise level is near the critical point to lose the robust evolutionary process. For 
a = 0.005, the decrease stops after 20 generations, and the variance values scatter 
for later generations. 




and so forth. On the other hand, the landscape evolved at a < ac is rugged. 
Except for the vicinity of the given initial conditions, the expression dynamics 
do not reach the target pattern. 

The observed proportionality between Vip and Vg is not self-evident. In- 
deed, if a random network is considered for a gene regulatory network, such 
proportionality is not observed. In the present simulation, after a few gen- 
erations of evolution, both the variances decrease, following proportionality, 
if (7 > (Tc- Although there is no complete derivation for this relationship, it 
is suggested that this proportionality as well as the relationship Vip > Vg 
is a consequence of evolutionary stability to keep a single-peakedness in the 
distribution P{x = phenotype,a = genotype), under conditions of strong 
selective pressure and low mutation rate jl5u35) . 



4.2 Proportionality between the two variances across genes 



As mentioned above, the gene expression dynamics evolved under a sufficient 
level of noise have a characteristic property; the attractor providing the phe- 



11 



notype of the highest fitness has a large basin volume and is, hence, attracted 
globally by a developmental process under noise. 

Note that the expression level Xj of non-target genes j could be either on 
or off, because there is no selection pressure directed at fixing their expression 
level. Still, each expression level Xj can have some correlation with the fitness 
in general. Hence it is also interesting to study the variance of each expression 
level and discuss its evolutionary changes. 

Similar to the variances for the fitness, the phenotypic variance Vip{i) for 
each gene i in an isogenic population is defined on the basis of the variance 
of the expression of each gene i, with each Xi — Sign{xi — di), in an isogenic 
population. Accordingly, the variance computed by using the distribution of 
Xi in this heterogenic population gives Vg{i) for each gene 

Following Price equation [21], the rate of change in each expression level 
between generations is expected to be correlated with Vg(i), as it has direct 
or indirect influence to the fitness, through the gene expression dynamics (1). 
In contrast, here, we are interested in the variance of isogenic phenotypic fluc- 
tuations of each expression level Vip{i). Indeed, this variance decreases over 
generations for most genes. As the evolution progresses, these expressions 
also start to be rigidly fixed so that their variances decrease over most genes. 
In Fig. 3(a), we have plotted Vg{i) versus Vip{i) over generations for several 
genes i. We can see that they decrease (roughly) in proportion, over gen- 
erations. Furthermore, the proportion coefficient Vg(i)/Vip{i) seems to take 
close values across many genes. 

In Fig. 3(b) we have plotted {Vip{i), Vg{i)), across all expressed genes, 
after evolution reached the genotype with the highest fitness. As shown, 
the proportionality (or strong correlation) between Vg{i) and Vip{i) holds 
across many (expressed) genes for a system through evolution. The ratio p = 
Vg{i) /Vip{i) increases with the decrease in cr, and at around ctc, it approaches 
- 1. 

Although the origin of the proportionality has not yet been completely un- 
derstood, a heuristic argument is proposed by using the distribution P{xi, a) 
for each gene expression Xi and by further assuming that the distribution 
maintains a single peak up to a common mutation rate (i.e., a common mu- 
tation rate for the error catastrophe threshold) over genes [ID]. The latter 
hypothesis may be justified as a highly robust system: In such a system with 
increased error-threshold mutation rates, once the error occurs, it propagates 
and percolates to many genes, so that a common error threshold value is ex- 
pected. Note that there are some preliminary experimental supports on this 
proportionality of the two variances over genes or phenotypic traits [T8ll46l 
1471148] . although future studies are required for the confirmation. 



^ Throughout the present paper, Vg{i) and Vip{i) with (i) denote the variances 
of each phenotype, expression level i, while, Vg and Vip denote the variances of the 
fitness. 
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Fig. 3 (a) Plot of {V,p{i),Vg{i)) for the genes i = 16,21,22,31,52 for the gener- 
ations 5-40. Each of the variances decreases over generations, so that variance for 
each gene changes from right (upper) to left (lower) in the figure. As described 
in the text, Vip{i) was computed as the variance of the distribution of Sign{xi) 
over L — 500 runs for an identical genotype, while Vg{i) was computed as a vari- 
ance of the distribution of {Sign{xi)) over = 500 individuals, where Sign{xi) 
refers to the mean over 500 runs. With generations, both the variances decrease 
roughly with proportion, with a trend of common proportion coefficient, (b) Plot 

of {Vip{i),Vg{i)) across all expressed genes i (i.e. for such genes i that x{i) > 9i), 
after evolution is completed (for the generations 25-60), for three values of noise 
levels a = 0.005 (*), 0.03 (x), and 0.06 {+). 
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5 Phenotypic Evolution under Environmental Variation 

5.1 Restoration of plasticity with the increase in fluctuations by 
environmental change 

So far, through the selection process under a fixed fitness/environmental con- 
dition, both the fluctuations and the rate of evolution decrease. The system 
loses plasticity against the change caused by external noise or external mu- 
tation. Nevertheless, in nature, neither the fluctuations nor the evolution 
potential vanish. How are phenotypic plasticity, fluctuations, and evolution- 
ary potential sustained in nature? 

One possible origin for the preservation of plasticity may be environmen- 
tal fluctuation [51], as has also been studied in terms of statistical Dhvsics[52l 
[5?1[53] . The plasticity of a biological system is relevant for coping with the en- 
vironmental change that may alter phenotypic dynamics in order to achieve 
a higher fitness. In the present modeling, there can be two ways to include 
such an environmental change. 

One is a direct method, in which an input term given in eq.(4) is changed 
with each generation, while the fitness condition is maintained. The other 
method is indirect, in which environmental change is introduced as the change 
in the fitness condition, while preserving the dynamics itself. Here, we discuss 
the simulation result of the former procedure first and will consider the result 
from the other procedure later in §5.3. 

To change the environmental condition, we varied the input pattern at 
some generation. Here, we change the input pattern Ij (1 < j < fcinp) from 
the generation at which the system had already adapted to the environment 
and decreased the phenotypic variances. An example is plotted in Fig. 4. Here, 
Ij initially takes Ij = 1 for 1 < j < kinp/2 and Ij — for kinp/2 < j < kmp, 
before switching to 1 — Ij after the 29th generation. By switching the envi- 
ronment, the fitness first decreases and later adapts to the new environment 
(Fig.4a). 

To determine the evolution of phenotypic plasticity, we computed the 
variances of the fitness, Vip and Vg, over successive generations (see Fig. 4b). 
After the switch of the fitness condition, both Vip and Vg first increase to a 
relatively high level and continuously increasing further over a few genera- 
tions. At later generations, both Vip and Vg again decrease, maintaining the 
proportionality. The proportionality law between the genetic and epigenetic 
variances is satisfied with both increase and decrease in plasticity through 
the evolution. 

With the increase in Vip, the fitness is more variable with noise, which also 
leads to higher changeability against environmental conditions. The gene ex- 
pression dynamics regain plasticity, which allows for the switch of the target 
genes after further generations. Then, with the increase in Vg, the change- 
ability against genetic change increases, thus increasing the evolvability. 

Next, we explore the change in the variances of the expression of each 
gene. As shown in Fig. 5, both variances Vip{i) and Vg{i) increase, as a result 
of environmental change. Expressions of all genes are more variable against 
noise and mutation. Most gene expressions gain higher plasticity by increas- 
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Fig. 4 The time course in the fitness and the variance of the fitness over genera- 
tions. First, the evolution under an environment Ij — 1 for 1 < J < 4, and Ij = 
for 5 < J < kinp = 8 is simulated to progress up to 30 generations. After the 29th 
generation, we switch the fitness condition to Ij = for 1 < j < 4, and Ij = 1 for 
5 < j < 8, which is maintained at later generations. M = 64, L = N = 700, and 6i 
is distributed uniformly in [0.1, 0.3]. The switch initially causes a decrease in the 
fitness, but after a few dozens of generations, almost all networks evolve to adapt 
to the new fitness condition. The noise level is set at ct = 0.06 > CTc- Top: The 
time course of the average fitness and the variance Vip throughout the evolution. 
Bottom: The plot of the variances of the fitness, Vg versus Vip for each generation 
after the switch of the environmental condition. The generations (up to 60) are 
represented by different colors. Both the variances increase in correlation, after the 
switch, and later, they decrease in proportion, to adapt to the new condition. 



15 




Fig. 5 Change in the variances Vip{i) and Vg{i) /Vip{i) at after the switch in en- 
vironmental condition after the 29th generation, as given in Fig. 4. The variance 
Vip{i)) is plotted up to generation 60, where the switch is given at the 30th genera- 
tion. Plotted for the genes i =13,15,28,40,47,51. For the genes 13, 28, and 47, Vg{i) 
before the switch after 29th generation is smaller than 10~^. The color represents 
Vg{i)/Vip{i), as shown in the right bar. 



ing the variance in their expression. Besides the increase in the variances, 
the ratio pi — Vg{i) /Vip{i) also increase at some generation, to approach 
Vip ~ Vg (see the color change in Fig. 5, where the red color shovifs higher 
ratio Pi. This increase in Vg{i)/Vip{i) suggests that the system is closer to 
the error catastrophe point, where the stability condition in the distribution 
function P{xi = phenotype, a = genotype) is lost. This leads to the increase 
in plasticity, with which the adaptation to a novel condition is achieved. 
Once the fitted phenotypes are generated by this adaptation, the variances 
decrease, with restoring the proportionality between Vg{i) and Vip{i). 

To sum up, adaptation to novel environment is characterized by the phe- 
notypic variances as follows. With the increase in Vip{i), the sensitivity of the 
phenotype to noise and environmental change is increased, thereby increasing 
plasticity. With the increase in Vg{i), changeability of the phenotype by mu- 
tation is increased, thereby accelerating evolutionary change of each expres- 
sion level. With the increase in Vg{i) /Vip{i), the sensitivity to genetic change 
is further increased, thus facilitating the evolution. With these trends — the 
increase in Vip{i), Vg{i), pi = Vg{i)/Vip{i), the adaptation to a novel environ- 
ment is fostered. With this increase in plasticity, gene expression dynamics 
for adapting to a novel condition are explored. Once these fitted dynamics 
are shaped, the variances decrease, leading to a decrease in plasticity and 
increase in robustness. 
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Fig. 6 The average of the mean fitness < F > plotted for each generation, under 
continual environmental variation, as described in the text. M = 64, k = 8, and 
kinp — 4, while Ij are changed randomly within [0,.8]. The average of the mean 
fitness, F, of each individual (over L = 100 runs) is computed over the total 
population (A'^ = 100) at each generation. The noise level, a, is 0.1 (red), 0.05 
(~ (Tc; green), O.Ol(blue), and 0.001 (pink). At around a ~ 0.05, the average fitness 
reaches the highest level. 



5.2 Optimal noise level for varying environment 

Now, we consider evolution under continual environmental variation. To dis- 
cuss a long-term environmental change, we switch the environmental condi- 
tion by generation. To be specific, we change randomly li within [0, 1] per 
generation. 

When environmental changes are continuously repeated, the decrease and 
increase in the variances Vip and Vg are repeated. Note that it takes more 
generations to adapt to a new fitness condition, if the phenotypic variances 
have been smaller. In our model, if the noise level in development is larger, 
the phenotypic variances already take a small value during the adaptation 
to satisfy the fitness condition. Hence, in this case, it takes more generations 
to adapt to a new fitness condition. On the other hand, if the noise level 
cr is smaller than CTc ^ -05, robust evolution does not progress. Hence, for 
continuous environmental change, there will be an optimal noise level to 
both adapt sufficiently fast to a new environment and evolve the robustness 
of fitness for each environmental condition. In Fig. 6, we have plotted the 
time course of the average fitness in population. If the noise level is large, 
the system cannot follow the frequent environmental change and the average 
fitness cannot increase sufficiently. On the other hand, if the noise level is 
small, the fitness increases; however, if it is too small, the fitness of some 
individuals remains rather low. Indeed, there is an optimal noise level at 
which the average fitness is maximal, as shown in Fig. 7, where the average 
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Fig. 7 The average fitness and the variances, Vip and Vg, through the course of the 
evolution under environmental variation as in Fig. 6. Each value is further averaged 
over 200-700 generations. The overall temporal averages are plotted against the 
noise level a. Instead of the fitness itself, its sign inversion — < _F > is plotted. At 
around a ~ 0.05, the average fitness takes a maximal value, and at a slightly below 
it, Vg approaches Vip. 
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Fig. 8 The variances of the fitness {Vip, Vg) through the course of the evolution 
under environmental variation as in Fig. 6. Each point is a result of one generation, 
and the plot is taken over 200-700 generations. The noise level ct is 0.1 (red), 0.005 
(~ (Tc; green), 0.01 (blue), and O.OOl(pink). 
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fitness over generations is plotted against the noise level a. This optimal noise 
level is close to the value of the robustness transition ac- 

Next, we plotted the variances Vip and Vg over generations (see Fig. 8). 
When a < ac, then Vg > Vip and both the variances remain rather large, 
demonstrating that robustness has not evolved at all. For a ac, Vg < Vip 
and the variances remain small. The robustness has evolved, but the system 
cannot adapt to an environmental change as the variances have become too 
small. In contrast, for a ^ ac, Vip and Vg vary between low and high values 
over generations, maintaining the proportionality between the two variances, 
with Vip slightly larger than Vg. 



5.3 Adaptation against switches of the fitness condition 
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Fig. 9 The average of fitness plotted per generation, where the fitness condition 
in the model (l)-(2) is switched for every 10 generations between + + + + + + ++ 

and + + + H . The average of the mean fitness F of each individual (over 

L = 200 runs) is computed over the total population (A^ — 200) at each generation. 
The noise level a is 0.1 (red), 0.008 (~ a^; green) and 0.001 (blue). 

For confirmation of the result in the last section, we also carried out 
numerical experiments by adopting a separate procedure, i.e., by switching 
the fitness condition. As a specific example, we carried out the simulation by 
taking the model (l)-(2), without including input terms (4). After the gene 
expression dynamics are evolved with the fitness to prefer > for the 
target genes i = 1,2, ...k{= 8) adopted already, then at a certain generation, 
we change the fitness condition so that the genes i — 1,2, .., k/2 are on and 

the rest are off (i.e., the fittest gene expression pattern is + + + H , 

instead of + + + + + + ++: In this model, the gene is off if a;^ < 0). Here, 
we switch after sufficiently large generations when the fittest networks are 
evolved (i.e., with Xi > for target genes). 
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Fig. 10 The temporal average of the average fitness and the variances Vip and Vg 
in the evolution, with the change in the fitness condition for every 10 generations 
as in Fig. 9. The overall temporal averages over population and over generations 
are plotted against the noise level a. Instead of the fitness itself, its sign inversion 
— < F > is plotted. The temporal average is taken over 500-1000 generations. At 
around a ~ .008, the average fitness takes a maximal value, and at a slightly below 
it, Vg exceeds Vip. 
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Fig. 11 Variances of the fitness, (Vp, Vg), are plotted over generations through the 
course of the evolution with fitness change after every 10 generations, as described 
in Fig. 9. Each point is a result of one generation, and the plot is taken over 500- 
1000 generations. The noise level a is 0.1 (red), 0.02(green), .008 (~ (Tc;blue), and 
0.001 (pink). 



In this case as well, the variances of the fitness, Vi^ and Vg, first increase 
in proportion to adapt to a new fitness condition. Later, they decrease in pro- 
portion to gain robustness to noise and mutation. Next, we again computed 
Vip(i) and Vg(i) successively through the course of the evolution. Immedi- 
ately after the switch in the fitness, the variances Vip(i) and lg(i) increase 
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as well as the ratio Vg{i)/Vip{i). These variances in gene expression levels 
facilitate plasticity, and adaptation to a new environment. 

When environmental changes are continuously repeated, the decrease and 
increase processes of the variances Vip and Vg are repeated. In Fig. 9, we have 
plotted the time course of the average fitness in population, when the fitness 
condition is switched after every 10 generations. In this case again, when 
the noise level is near Uc, fast adaptation to a new environment and the 
increase of robustness at later generations are compatible. Dependence of 
the average fitness on the noise level is shown in Fig. 10, which also shows 
an optimal noise level near the; robustness transition ctc (which is lower than 
the case in §5.2). Indeed, below this noise level, Vg exceeds Vip and the 
robustness to mutation is lost. The plot of the variances Vip and Vg over 
generations in Fig. 11 shows that at cr ^ CTc, they go up and down, maintaining 
an approximate proportionality between the two, with Vip slightly less than 
Vg. Overall, the behavior here under the fitness switch agrees well with that 
under the environmental variation in §5.2. 

6 Summary and Discussion 

In the present paper, we have studied biological robustness and plasticity, 
in terms of phenotypic fluctuations. The results are summarized into three 
points. 

(1) Confirmation of earlier results on the phenotypic variances due to 
noise and due to genetic variation: The two variances, Vip (due to noise) 
and Vg (due to mutation) decrease in proportion through the course of ro- 
bust evolution under a fixed environmental condition. After the evolution to 
achieve robustness to noise and mutation is completed, Vip{i) and Vg{i) are 
proportional across expressions of most genes (or different phenotypic traits). 
In short, 

plasticity (changeability) of phenotype oc V^p oc oc evolution 
speed 

through the course of the evolution and across phenotypic traits (expressions 
of genes). 

(2) Increase in the phenotypic variances and recovery of plasticity: When 
robustness is increased under a given environmental condition, the system 
loses plasticity to adapt to a novel environment. When the environmental 
condition is switched, both the phenotypic variances due to noise and due to 
genetic variation increase to gain plasticity and thus to adapt to the novel 
environment. This increase is observed both for the variance of fitness and 
of each gene expression level. 

(3) Optimal noise level to achieve both robustness and plasticity under 
continuous environmental change. There is generally a threshold level of noise 
in gene expression (or in developmental dynamics), beyond which robustness 
to noise and mutation evolves. If the noise level is larger, however, the system 
loses plasticity to adapt to environmental changes, whereas if it is much lower, 
a robust, fitted phenotype is not generated. At around the noise level for the 
" robustness transition," the system can adapt to environmental changes and 
achieve a higher fitness. There, the phenotypic variances Vip and Vg increase 
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and decrease, roughly maintaining the proportionahty between the two, while 
sustaining Vip ^ Vg. 

From a statistical-physicist viewpoint, the relationship between respon- 
siveness and fluctuation is expected as a proper extension of the fluctuation- 
response relationship. In this context, it is rather natural that the response ra- 
tio to environmental change, i.e., plasticity, is proportional to Vip, the isogenic 
phenotypic fluctuation. On the other hand, as Fisher stated, the variance 
due to genetic change, Vg, is proportional to evolution speed, i.e., response 
of the fitness against mutation and selection. Interestingly, our simulations 
and evolutionary stability argument suggest the proportionality between Vip 
and Vg. This implies the proportionality between environmental plasticity 
and evolvability. 

In fact, Waddington |49[[5n] coined the term genetic assimilation, in which 
phenotypic changes induced by environmental changes foster later genetic 
evolution. Since then, positive roles of phenotypic plasticity in evolution have 
been extensively discussed [31[11I2]. Our study gives a quantitative representa- 
tion of such relationship in terms of fluctuations. 

Existence of the threshold noise level below which robustness is lost is 
reminiscent of a glass transition in physics; For a higher noise level, dynam- 
ical systems for global attraction to a functional phenotype are generated 
through evolution, whereas for a lower noise level, the dynamics follow mo- 
tion in a rugged landscape, where perturbation to it leads to a failure in 
the shaping of the functional phenotype. In fact, Sakata et al. considered a 
spin-glass model whose interaction matrix evolves to generate a high-fitness 
thermodynamic state. The transition to lose robustness was found by lower- 
ing the temperature. Interestingly, this transition is identified as the replica 
symmetry breaking transition from a replica symmetric phase 33,34 . 

Under environmental fiuctuation, the evolution to achieve both plasticity 
to a new environment and robustness of a fitted state is possible near this 
transition, for losing the robustness. In other words, one may regard that a 
biological systems favors the "edge-of- glass" state. 
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